O exercício consiste em construir um ou mais modelos lineares em uma base pública de dados sobre nascimentos na Carolina do Norte (EUA) no ano de 2001, que consiste em uma amostra de 1450 registros de nascimento selecionados pelo estatístico John Holcomb.
Carregue a base de dados NCBirths, disponível em https://vincentarelbundock.github.io/Rdatasets/datasets.html.
Desenvolva os itens a seguir aplicando técnicas gráficas e formais, e apresente resultados, explicações e considerações que julgar necessários.
Utilize o presente documento Markdown (referências em http://rmarkdown.rstudio.com, https://bookdown.org/yihui/rmarkdown/) para a inserção do código e geração de resultados.
Regressão Múltipla
# bloco de código - item a
NCbirths = read.csv("NCbirths.csv")
nc_births = NCbirths
str(nc_births)## 'data.frame': 1450 obs. of 16 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Plural : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Sex : int 1 2 1 1 1 1 2 2 2 2 ...
## $ MomAge : int 32 32 27 27 25 28 25 15 21 27 ...
## $ Weeks : int 40 37 39 39 39 43 39 42 39 40 ...
## $ Marital : int 1 1 1 1 1 1 1 2 1 2 ...
## $ RaceMom : int 1 1 1 1 1 1 1 1 1 1 ...
## $ HispMom : Factor w/ 6 levels "C","M","N","O",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Gained : int 38 34 12 15 32 32 75 25 28 37 ...
## $ Smoke : int 0 0 0 0 0 0 0 0 0 0 ...
## $ BirthWeightOz: int 111 116 138 136 121 117 143 113 120 124 ...
## $ BirthWeightGm: num 3147 3289 3912 3856 3430 ...
## $ Low : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Premie : int 0 0 0 0 0 0 0 0 0 0 ...
## $ MomRace : Factor w/ 4 levels "black","hispanic",..: 4 4 4 4 4 4 4 4 4 4 ...
summary(nc_births)## X ID Plural Sex
## Min. : 1.0 Min. : 1.0 Min. :1.000 Min. :1.000
## 1st Qu.: 363.2 1st Qu.: 363.2 1st Qu.:1.000 1st Qu.:1.000
## Median : 725.5 Median : 725.5 Median :1.000 Median :1.000
## Mean : 725.5 Mean : 725.5 Mean :1.037 Mean :1.487
## 3rd Qu.:1087.8 3rd Qu.:1087.8 3rd Qu.:1.000 3rd Qu.:2.000
## Max. :1450.0 Max. :1450.0 Max. :3.000 Max. :2.000
##
## MomAge Weeks Marital RaceMom HispMom
## Min. :13.00 Min. :22.00 Min. :1.000 Min. :1.000 C: 2
## 1st Qu.:22.00 1st Qu.:38.00 1st Qu.:1.000 1st Qu.:1.000 M: 128
## Median :26.00 Median :39.00 Median :1.000 Median :1.000 N:1283
## Mean :26.76 Mean :38.62 Mean :1.345 Mean :1.831 O: 3
## 3rd Qu.:31.00 3rd Qu.:40.00 3rd Qu.:2.000 3rd Qu.:2.000 P: 9
## Max. :43.00 Max. :45.00 Max. :2.000 Max. :8.000 S: 25
## NA's :1
## Gained Smoke BirthWeightOz BirthWeightGm
## Min. : 0.0 Min. :0.0000 Min. : 12.0 Min. : 340.2
## 1st Qu.:20.0 1st Qu.:0.0000 1st Qu.:106.0 1st Qu.:3005.1
## Median :30.0 Median :0.0000 Median :118.0 Median :3345.3
## Mean :30.6 Mean :0.1446 Mean :116.2 Mean :3295.6
## 3rd Qu.:40.0 3rd Qu.:0.0000 3rd Qu.:130.0 3rd Qu.:3685.5
## Max. :95.0 Max. :1.0000 Max. :181.0 Max. :5131.4
## NA's :40 NA's :5
## Low Premie MomRace
## Min. :0.00000 Min. :0.0000 black :332
## 1st Qu.:0.00000 1st Qu.:0.0000 hispanic:164
## Median :0.00000 Median :0.0000 other : 48
## Mean :0.08621 Mean :0.1317 white :906
## 3rd Qu.:0.00000 3rd Qu.:0.0000
## Max. :1.00000 Max. :1.0000
##
Uma inspeção visual já permite identificar NAs nos campos Weeks, Gained e Smoke, mas os campos que a descrição diz serem categóricos foram tratados como numéricos, e seus sumários não fazem sentido.
Não foi especificado um problema para ser estudado. Vejamos como cada campo poderia ser interpretado se fosse a variável dependente. Há dois grupos: Características da mãe e as da gestação e bebê:
MomAge: Embora se apresente como numérica discreta, pode ser considerada contínua.Marital: Categórica com apenas dois valores (Casada ou não)MomRace, RaceMom e HispMom, todos categóricosSmoke, categóricoGained, apesar de ser o peso ganho pela mãe, foi considerado parte da gestação e comentado na seção seguintePlural: Categórica, pode ser uma logística por Single/Not singleSex: Categórica, também pode ser uma logística por Male/FemaleGained: ContínuaWeeks: Numérica discretaPremie: Categórica, depende de WeeksBirthWeightOz ou BirthWeightGm, numéricas. Em Oz, parece ser apenas inteiro enquanto em g, é float.
Plural!= “Single”)Low: Categórica. Adequada para logística. Depende integralmente do peso do bebê.MomRace, HispMom e MomRace (ô povo racista!): Há grande sobreposição entre os conceitos.
Weeks e Premie: A segunda depende completamente da primeiraBirthWeight{Oz,Gm} e Low: Mesmo caso: A segunda também depende completamente das primeiras
Low realmente têm peso menor que 2500gTodas as colunas poderiam ser sujeitas a tentativas de explicação, mas suas características são diferentes:
| Quantitativas | Categóricas | |
|---|---|---|
| Mãe | Gained¹, MomAge |
Smoke, Marital, “raças”³ |
| Gestação | Gained¹, Weeks |
Plural, Premie |
| Bebê | Peso² | Sex, Low |
As “raças” estão espalhadas por três variáveis: MomRace, RaceMom e HispMom
As variáveis quantitativas podem ser objeto de estudo de modelos lineares, enquanto as categóricas precisam ser explicadas por outros modelos. Em particular, as binárias (Sex, Smoke, Marital, Premie e Low) poderiam ser estudadas com a regressão logística.
Num contexto humano, não faz sentido prever grupo étnico, estado civil ou idade da mãe com base nos dados disponíveis, então as variáveis MomAge, Marital e as “raças” não serão objeto de estudo dessa natureza.
1 - Quem ganha o peso é a mãe, mas isso decorre da gravidez
2 - O peso é dado em duas unidades, basta ignorar uma delas
3 - Geneticamente, o conceito de raça humana não existe, mas culturalmente a distinção pode ser importante, como no caso do local da origem desses dados.
As três variáveis quantitativas que podem (talvez) serem estudadas com regressão linear múltipla são:
Gained)BirthWeightOz)Weeks)hist(nc_births$Gained, density = 10, breaks = 30)# density(nc_births$Gained)
sw_gained = shapiro.test(nc_births$Gained)
gained_mean = mean(nc_births$Gained, na.rm = T)A distribuição das frequências de pesos ganhos parece ser normal, com média 30.6014184 e desvio padrão 13.8774929 (apenas os valores presentes. 40 estão faltando), o resultado de um teste Shapiro-Wil é 0.9843275 com \(p < 5\%\) (3.12516610^{-11}).
TODO: adasad
peso_normal = nc_births[nc_births$Low==0,]
baixo_peso = nc_births[nc_births$Low==1,]
par(mfrow = c(1, 2))
boxplot(peso_normal$BirthWeightOz,
main="Bebês com peso acima de 2500g")
boxplot(baixo_peso$BirthWeightOz,
main="Bebês com peso baixo")par(mfrow = c(1, 1))
hist(nc_births$BirthWeightOz, breaks=25,
# xlab="Semanas de gestação",
main="Peso do bebê (geral)")hist(peso_normal$BirthWeightOz, breaks=25,
main="Peso normal")hist(baixo_peso$BirthWeightOz, breaks=25,
main="Baixo Peso")A distribuiçã das frequências de duração da gravidez é enviesada em direção à duração esperada da gestação a termo: 40 semanas após a última menstruação, 38 após a fertilização, e não pode ser considerada uma distribuição Normal (Teste Shapiro-Wilk: W=0.8506, \(p=6.3253773\times 10^{-35}\)).
A média é 38.621 e o desvio padrão é de 2.699.
O peso do bebÊ é dado tanto em onças (Oz) quanto em gramas (g). Se não houver erros de conversão de unidade de massa, basta eliminar uma delas, pois têm rigorosamente a mesma informação.
Uma eventual divergência, no entanto, pode indicar erros graves na coleta e/ou registro dos dados, e uma escolha mais complexa precisa ser feita.
# bloco de código - item b
summary(nc_births$BirthWeightGm/nc_births$BirthWeightOz)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 28.35 28.35 28.35 28.35 28.35 28.35
if (mean(nc_births$BirthWeightGm/nc_births$BirthWeightOz)/max(nc_births$BirthWeightGm/nc_births$BirthWeightOz) != 1 ) {
error("Há um erro de conversão entre Oz e g")
} else {
if (max(nc_births$BirthWeightGm/nc_births$BirthWeightOz) == 28.35) {
print("Todas as conversões de unidade estão corretas")
} else {
print("Embora consistentes, as conversões não estão corretas")
}
}## [1] "Todas as conversões de unidade estão corretas"
Conclui-se que as colunas BirthWeightGm e BirthWeightOz contêm informação idêntica, e uma delas pode ser ignorada com segurança.
Parece haver uma confusão no registro de raças. Há três campos para isso, e encontrei pelo menos uma inconsistência no aparente registro de mães hispânicas como japonesas:
# for (race in unique(nc_births$MomRace)) { print(paste(race, nrow(nc_births[nc_births$MomRace == race,])))} # TODO Um Summary não seria melhor?
summary(nc_births$MomRace)## black hispanic other white
## 332 164 48 906
summary(nc_births[nc_births$MomRace=="hispanic", c("RaceMom", "HispMom")]) # Mostra que `RaceMom` é 5 para todas com `MomRace == hispanic`## RaceMom HispMom
## Min. :5 C: 2
## 1st Qu.:5 M:128
## Median :5 N: 0
## Mean :5 O: 3
## 3rd Qu.:5 P: 8
## Max. :5 S: 23
JapasSoQueNao = nrow(nc_births[nc_births$RaceMom==5,]) # "Japanese", só que não
print(paste("Registros marcados como mães japonesas: ", JapasSoQueNao))## [1] "Registros marcados como mães japonesas: 164"
print("Mães com RaceMom == 5")## [1] "Mães com RaceMom == 5"
summary(nc_births[nc_births$RaceMom==5, c("MomRace", "HispMom")])## MomRace HispMom
## black : 0 C: 2
## hispanic:164 M:128
## other : 0 N: 0
## white : 0 O: 3
## P: 8
## S: 23
Aparentemente, os 164 registros de mães “japonesas” na verdade são hispânicas com a marca incorreta no campo RaceMom mas correto em MomRace. O que me leva a crer que o correto seria considerar todas como hispânicas é a variedade do campo HispMom, que tem apenas 3 marcadas como “Other hispanic” (e que ainda não sei se há alguma chance de serem japonesas), e a ausência de uma marca “Not hispanic”.
Low está correto?O campo Low, segundo a descrição, chama a atenção para bebês cujo peso ao nascer é inferior a 2500g.
# summary(nc_births[nc_births$Low == "Y", c("BirthWeightOz", "BirthWeightGm")])
# summary(nc_births[nc_births$Low == "N", c("BirthWeightOz", "BirthWeightGm")])
# summary(nc_births[nc_births$BirthWeightGm > 2500, c("Low", "BirthWeightOz", "BirthWeightGm")])
nrow(nc_births[nc_births$BirthWeightGm > 2500 & nc_births$Low=='Y', c("Low", "BirthWeightOz", "BirthWeightGm")]) # Deve ser zero## [1] 0
nrow(nc_births[nc_births$BirthWeightGm <= 2500 & nc_births$Low=='N', c("Low", "BirthWeightOz", "BirthWeightGm")]) # Deve ser zero## [1] 0
Foi confirmado que, realmente, os registros marcados como Low==Y todos têm peso abaixo de 2500g, e os marcados com Low==N todos têm peso acima do informado na documentação. De forma simétrica, nenhum registro inconsistente (isto é, com marca incorreta) foi encontrado, o que caracteriza uma dependência completa.
Por isso, esse campo pode seguramente ser ignorado nas regressões.
Por inspeção visual, foi possível detectar valores ausentes nas colunas Weeks(1), Smoke(5) e Gained(40). Faremos avaliação de cada caso assim que der.
col.sums = colSums(is.na(nc_births))
col.sums[col.sums>0] # Verificação programática de quais colunas têm NAs## Weeks Gained Smoke
## 1 40 5
nc_births<- cbind(nc_births, somaNulos = rowSums(is.na(nc_births)))
NCB_NAs = nc_births[nc_births$somaNulos > 0,]
NCBnoNAs = nc_births[complete.cases(nc_births),]Parece preocupante a grande quantidade de ausências em Gained, que corresponde a 2.7586207% das linhas. A estratégia ideal depende da alavancagem desses registros.
Segundo a documentação do DataSet, várias dessas colunas podem ser convertidas em fatores:
# Fatorizar o que der
source("factorize.R")
useEtn = TRUE # Flag para unificar campos étnicos
nc_births = NCB_noId(nc_births)
nc_births = NCB_factorize(nc_births, useEtn)
ncb_births_noGainedNA = nc_births[!is.na(nc_births$Gained), ]
ncb_births_noGainedNA = ncb_births_noGainedNA[ncb_births_noGainedNA$somaNulos == 0, ]
ncb_births_noNAsAtAll = nc_births[nc_births$somaNulos == 0, ]
if (useEtn) {
COLS_W = c("Plural", "Sex", "MomAge", "Weeks", "Marital",
"Gained", "Smoke", "BirthWeightOz", "Premie", "Etnicidade")
} else {
COLS_W = c("Plural", "Sex", "MomAge", "Weeks", "Marital",
"RaceMom", "HispMom", "Gained", "Smoke",
"BirthWeightOz", "Premie", "MomRace")
}
ncb_weight = nc_births[, COLS_W]
# ncb_weight_noNAs = ncb_births_noGainedNA[, COLS_W] # ncb_births_noNAsAtAll
ncb_weight_noNAs = ncb_births_noNAsAtAll[, COLS_W] #
ncb_weigh_patchedNA = ncb_weight
# ncb_weigh_patchedNA[ncb_weigh_patchedNA$somaNulos >0, "Gained"] = mean(ncb_weigh_patchedNA$Gained, na.rm = T)
ncb_weigh_patchedNA[is.na(ncb_weigh_patchedNA$Gained), "Gained"] = mean(ncb_weigh_patchedNA$Gained, na.rm = T)
ncb_weigh_patchedNA[is.na(ncb_weigh_patchedNA$Smoke), "Smoke"] = "N"
# $Gained[is.na(ncb_weight$Gained)] = mean(ncb_weight$Gained)
# Remover coluna desnecessária
ncb_births_noGainedNA$somaNulos <- NULL
nc_births$somaNulos <- NULL# boxplot(nc_births$Gained ~ nc_births$MomRace + nc_births$RaceMom, main='Peso ganho por "Raça"') # Pouco útil sem fatoresalias(lm(BirthWeightOz ~ . - BirthWeightGm, data = new_ncbirths(NCbirths)))## Model :
## BirthWeightOz ~ (Plural + Sex + MomAge + Weeks + Marital + RaceMom +
## HispMom + Gained + Smoke + BirthWeightGm + Low + Premie +
## MomRace) - BirthWeightGm
##
## Complete :
## (Intercept) PluralTwins PluralTriplets SexF MomAge Weeks
## MomRacehispanic 0 0 0 0 0 0
## MomRaceother 0 0 0 0 0 0
## MomRacewhite 1 0 0 0 0 0
## MaritalNot Married RaceMomBlack RaceMomAmericanIndian
## MomRacehispanic 0 0 0
## MomRaceother 0 0 1
## MomRacewhite 0 -1 -1
## RaceMomChinese RaceMomHispanic RaceMomFilipino
## MomRacehispanic 0 1 0
## MomRaceother 1 0 1
## MomRacewhite -1 -1 -1
## RaceMomOtherAsianOrPacific HispMomM HispMomN HispMomO
## MomRacehispanic 0 0 0 0
## MomRaceother 1 0 0 0
## MomRacewhite -1 0 0 0
## HispMomP HispMomS Gained SmokeY LowY PremieY
## MomRacehispanic 0 0 0 0 0 0
## MomRaceother 0 0 0 0 0 0
## MomRacewhite 0 0 0 0 0 0
alias(lm(BirthWeightOz ~ . - BirthWeightGm, data = new_ncbirths(NCbirths, T)))## Model :
## BirthWeightOz ~ (Plural + Sex + MomAge + Weeks + Marital + Gained +
## Smoke + BirthWeightGm + Low + Premie + Etnicidade) - BirthWeightGm
A sobreposição de informações dadas nas colunas étnicas aparece no relatório criado pela função alias(). Isso pôde ser corrigido ao se unificar as informações em uma única coluna.
Curiosamente, há combinações inesperadas ali, como mães negras oriundas de Porto Rico e da América do Sul sendo consideradas também hispânicas. Por outro lado, distinguir índios de continentes distintos parece compreensível.
Tentativa de explicar o peso ganho em função das outras variáveis.
# summary(lm(Gained ~ ., data=nc_births)) # Ficou igual porque lm() eliminou as linhas sem valores para Gained
fit_gained1 = lm(Gained ~ ., data=ncb_births_noGainedNA)
summary(fit_gained1)##
## Call:
## lm(formula = Gained ~ ., data = ncb_births_noGainedNA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.624 -8.566 -0.883 7.691 61.092
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.88006 8.69879 2.285 0.022441 *
## PluralTwins 12.59628 2.26949 5.550 3.41e-08 ***
## PluralTriplets 24.00970 6.80442 3.529 0.000432 ***
## SexF 0.21819 0.71576 0.305 0.760542
## MomAge -0.14900 0.06712 -2.220 0.026589 *
## Weeks -0.10314 0.20964 -0.492 0.622795
## MaritalNot Married 0.79885 0.92547 0.863 0.388188
## SmokeY 0.87423 1.06189 0.823 0.410491
## BirthWeightOz 0.17957 0.02407 7.461 1.50e-13 ***
## BirthWeightGm NA NA NA NA
## LowY -0.78395 1.81751 -0.431 0.666293
## PremieY 1.75629 1.55438 1.130 0.258716
## Etnicidadeblack Black N -4.60985 3.00060 -1.536 0.124691
## Etnicidadeblack portoriq -18.38166 13.66639 -1.345 0.178837
## Etnicidadeblack south american -12.87652 13.62680 -0.945 0.344854
## Etnicidadecentro-south american 0.76697 4.07141 0.188 0.850607
## Etnicidadechinese 10.92400 9.86176 1.108 0.268178
## Etnicidadecuban -12.99595 9.84362 -1.320 0.186973
## Etnicidadefilipino 5.12258 13.65153 0.375 0.707540
## Etnicidadehispanic other -6.97672 8.23266 -0.847 0.396894
## Etnicidademexican -9.24698 3.16379 -2.923 0.003526 **
## EtnicidadeOtherAsianOrPacific -1.00753 4.08095 -0.247 0.805034
## Etnicidadeportoriq -0.22612 5.83565 -0.039 0.969096
## Etnicidadesouth americanIndian 19.21733 13.62796 1.410 0.158723
## Etnicidadewhite -2.36231 2.96472 -0.797 0.425699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.29 on 1385 degrees of freedom
## Multiple R-squared: 0.0976, Adjusted R-squared: 0.08261
## F-statistic: 6.513 on 23 and 1385 DF, p-value: < 2.2e-16
Curiosamente, vários campos que eu marquei para exclusão aparecem como NA no sumário, mas não todos: Low foi considerado, mas ficou terrivelmente insignificante.
Esse modelo explica apenas 8% do resultado. Será que existem modelos melhores?
# pairs(Gained ~ ., data=nc_births)
pairs(Gained ~ Plural + MomAge + Weeks + BirthWeightOz,
data=nc_births,
main = "Correlações entre os campos mais significativos"
)# fit_gained1_noNAs = lm(Gained ~ ., data=nc_births, na.action = na.omit)
# summary(fit_gained1_noNAs)
# NCBnoNAs = nc_births[]Não parece haver correlação nenhuma entre esses campos! O único padrão relevante é a já conhecida relação entre a idade gestacional do nascituro e o peso ao nascer: quanto mais próximo do termo (40 semanas), mais pesado; depois do termo, nem tanto.
Com base no sumário do modelo naïve, pode-se ver que os campos Plural e BirthWeightOz foram os únicos significativos.
fit_gained2 = lm(Gained ~ Plural + BirthWeightOz, data = ncb_births_noGainedNA)
summary(fit_gained2)##
## Call:
## lm(formula = Gained ~ Plural + BirthWeightOz, data = ncb_births_noGainedNA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.559 -8.976 -0.862 7.714 61.922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.71609 2.06272 5.680 1.64e-08 ***
## PluralTwins 12.46352 2.19662 5.674 1.69e-08 ***
## PluralTriplets 24.14434 6.79263 3.554 0.000391 ***
## BirthWeightOz 0.15823 0.01723 9.184 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.43 on 1405 degrees of freedom
## Multiple R-squared: 0.06504, Adjusted R-squared: 0.06304
## F-statistic: 32.58 on 3 and 1405 DF, p-value: < 2.2e-16
anova(fit_gained1, fit_gained2)## Analysis of Variance Table
##
## Model 1: Gained ~ Plural + Sex + MomAge + Weeks + Marital + Smoke + BirthWeightOz +
## BirthWeightGm + Low + Premie + Etnicidade
## Model 2: Gained ~ Plural + BirthWeightOz
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1385 244626
## 2 1405 253451 -20 -8825 2.4982 0.0002696 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Como resultado, todos os coeficientes são altamente significativos, e a análise de variância informa uma queda de mais de 8000 pontos na soma dos quadrados dos resíduos. Entretanto, o poder de explicação caiu para 6,34%.
Concluo que não vale a pena tentar prever o peso ganho durante a gestação com base nssses dados.
Da mesma forma que no modelo 1, vejamos com o que se parece um modelo com “tudo” dentro. Entretanto, como há campos com forte correlação entre si e com a variável dependente, precisamos nos livrar antes de:
Além disso, os 40 casos de ganho de peso ausentes precisam ser avaliados para uma tomada de decisão.
# ncb_weight = nc_births[,c("Plural", "Sex", "MomAge", "Weeks", "Marital",
# "RaceMom", "HispMom", "Gained", "Smoke",
# "BirthWeightOz", "Premie", "MomRace")]
#
# ncb_weight_noNAs = ncb_births_noGainedNA[c("Plural", "Sex", "MomAge", "Weeks", "Marital",
# "RaceMom", "HispMom", "Gained", "Smoke",
# "BirthWeightOz", "Premie", "MomRace"),]
# summary(lm(BirthWeightOz ~ ., data=ncb_weight, na.action = na.omit))
# fit_weight_1 é o fModelo 3 ====
fit_weight_1 = lm(BirthWeightOz ~ ., data=ncb_weight, na.action = na.omit)##
## Call:
## lm(formula = BirthWeightOz ~ ., data = ncb_weight, na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.588 -10.165 -0.473 10.194 51.977
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -31.76600 10.41147 -3.051 0.002324 **
## PluralTwins -24.91791 2.70283 -9.219 < 2e-16 ***
## PluralTriplets -32.97522 8.37525 -3.937 8.65e-05 ***
## SexF -3.29641 0.87689 -3.759 0.000178 ***
## MomAge 0.38212 0.08229 4.644 3.75e-06 ***
## Weeks 3.46722 0.24046 14.419 < 2e-16 ***
## MaritalNot Married -1.84540 1.14051 -1.618 0.105881
## Gained 0.28078 0.03227 8.701 < 2e-16 ***
## SmokeY -7.09110 1.29607 -5.471 5.30e-08 ***
## PremieY -7.48998 1.88547 -3.972 7.48e-05 ***
## Etnicidadeblack Black N -2.35540 3.70339 -0.636 0.524873
## Etnicidadeblack portoriq 7.33369 16.86260 0.435 0.663696
## Etnicidadeblack south american 4.18383 16.81073 0.249 0.803492
## Etnicidadecentro-south american 0.49574 5.02110 0.099 0.921366
## Etnicidadechinese 3.17904 12.16810 0.261 0.793931
## Etnicidadecuban 3.71059 12.14769 0.305 0.760064
## Etnicidadefilipino -30.03061 16.81461 -1.786 0.074320 .
## Etnicidadehispanic other 3.68649 10.15569 0.363 0.716662
## Etnicidademexican 2.36859 3.91357 0.605 0.545129
## EtnicidadeOtherAsianOrPacific -1.72621 5.03311 -0.343 0.731673
## Etnicidadeportoriq -13.26520 7.18070 -1.847 0.064911 .
## Etnicidadesouth americanIndian 7.88998 16.81736 0.469 0.639031
## Etnicidadewhite 1.41976 3.65717 0.388 0.697918
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.39 on 1386 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.4597, Adjusted R-squared: 0.4512
## F-statistic: 53.61 on 22 and 1386 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = BirthWeightOz ~ ., data = ncb_weigh_patchedNA, na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.250 -10.381 -0.569 10.246 52.373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.91916 10.24067 -3.605 0.000323 ***
## PluralTwins -23.90895 2.65762 -8.996 < 2e-16 ***
## PluralTriplets -32.33524 8.42338 -3.839 0.000129 ***
## SexF -3.41657 0.87099 -3.923 9.17e-05 ***
## MomAge 0.35694 0.08148 4.381 1.27e-05 ***
## Weeks 3.62095 0.23489 15.415 < 2e-16 ***
## MaritalNot Married -2.40859 1.12750 -2.136 0.032830 *
## Gained 0.28113 0.03244 8.666 < 2e-16 ***
## SmokeY -6.77780 1.29182 -5.247 1.78e-07 ***
## PremieY -6.62059 1.86564 -3.549 0.000400 ***
## Etnicidadeblack Black N -2.34155 3.72682 -0.628 0.529910
## Etnicidadeblack portoriq 7.00896 16.97642 0.413 0.679768
## Etnicidadeblack south american 4.47865 16.92608 0.265 0.791355
## Etnicidadecentro-south american 1.13483 5.00139 0.227 0.820532
## Etnicidadechinese 3.31510 12.25095 0.271 0.786738
## Etnicidadecuban 3.89491 12.23136 0.318 0.750201
## Etnicidadefilipino -29.98661 16.92982 -1.771 0.076736 .
## Etnicidadehispanic other 3.80635 10.22475 0.372 0.709748
## Etnicidademexican 2.49703 3.91547 0.638 0.523749
## EtnicidadeOtherAsianOrPacific -2.40474 5.01479 -0.480 0.631635
## Etnicidadeportoriq -9.51223 6.88705 -1.381 0.167441
## Etnicidadesouth americanIndian 7.62941 16.93276 0.451 0.652367
## Etnicidadewhite 1.11467 3.68020 0.303 0.762023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.5 on 1426 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4619, Adjusted R-squared: 0.4536
## F-statistic: 55.65 on 22 and 1426 DF, p-value: < 2.2e-16
As diferenças entre os coeficientes não são estatisticamente significantes:
comparaCoeficientes(sf1, sf1.1, "Coeficientes com e sem linhas incompletas")Chama a atenção o coeficiente para Etnicidadefilipino, com valor absoluto elevado e com significância de 0.055637 0.074320 (0.076736 no modelo que inclui as linhas com valores ausentes), perto do limite escolhido para \(\alpha\). Existe exatamente 1 registro de mãe com essa etnia.
if (useEtn) {
ncb_weight_NF = ncb_weight[ncb_weight$Etnicidade != "filipino",]
} else {
ncb_weight_NF = ncb_weight[ncb_weight$RaceMom != "Filipino",]
}
fit_weight_1b = lm(BirthWeightOz ~ ., data=ncb_weight_NF, na.action = na.omit)
sf1b = summary(fit_weight_1b)
sf1b##
## Call:
## lm(formula = BirthWeightOz ~ ., data = ncb_weight_NF, na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.588 -10.192 -0.478 10.198 51.977
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -31.76600 10.41147 -3.051 0.002324 **
## PluralTwins -24.91791 2.70283 -9.219 < 2e-16 ***
## PluralTriplets -32.97522 8.37525 -3.937 8.65e-05 ***
## SexF -3.29641 0.87689 -3.759 0.000178 ***
## MomAge 0.38212 0.08229 4.644 3.75e-06 ***
## Weeks 3.46722 0.24046 14.419 < 2e-16 ***
## MaritalNot Married -1.84540 1.14051 -1.618 0.105881
## Gained 0.28078 0.03227 8.701 < 2e-16 ***
## SmokeY -7.09110 1.29607 -5.471 5.30e-08 ***
## PremieY -7.48998 1.88547 -3.972 7.48e-05 ***
## Etnicidadeblack Black N -2.35540 3.70339 -0.636 0.524873
## Etnicidadeblack portoriq 7.33369 16.86260 0.435 0.663696
## Etnicidadeblack south american 4.18383 16.81073 0.249 0.803492
## Etnicidadecentro-south american 0.49574 5.02110 0.099 0.921366
## Etnicidadechinese 3.17904 12.16810 0.261 0.793931
## Etnicidadecuban 3.71059 12.14769 0.305 0.760064
## Etnicidadehispanic other 3.68649 10.15569 0.363 0.716662
## Etnicidademexican 2.36859 3.91357 0.605 0.545129
## EtnicidadeOtherAsianOrPacific -1.72621 5.03311 -0.343 0.731673
## Etnicidadeportoriq -13.26520 7.18070 -1.847 0.064911 .
## Etnicidadesouth americanIndian 7.88998 16.81736 0.469 0.639031
## Etnicidadewhite 1.41976 3.65717 0.388 0.697918
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.39 on 1386 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.4594, Adjusted R-squared: 0.4513
## F-statistic: 56.1 on 21 and 1386 DF, p-value: < 2.2e-16
comparaCoeficientes(sf1b, sf1)vif(fit_weight_1)## GVIF Df GVIF^(1/(2*Df))
## Plural 1.173561 2 1.040822
## Sex 1.007636 1 1.003811
## MomAge 1.315787 1 1.147078
## Weeks 2.134540 1 1.461006
## Marital 1.541343 1 1.241509
## Gained 1.050612 1 1.024994
## Smoke 1.099646 1 1.048640
## Premie 2.087255 1 1.444734
## Etnicidade 1.419620 13 1.013568
vif(fit_weight_1b)## GVIF Df GVIF^(1/(2*Df))
## Plural 1.173533 2 1.040816
## Sex 1.006888 1 1.003438
## MomAge 1.313160 1 1.145932
## Weeks 2.134147 1 1.460872
## Marital 1.540767 1 1.241276
## Gained 1.050590 1 1.024983
## Smoke 1.099513 1 1.048576
## Premie 2.087037 1 1.444658
## Etnicidade 1.415500 12 1.014584
As colunas Weeks e Premie são as que mais inflacionam a variância desse modelo. A retirada da mãe filipina diminuiu levemente a inflação do campo Etnicidade.
Menos é mais (de novo)
Baseado no contexto e nos valores observados das significâncias dos coeficientes do modelo 3, seleciono um conjunto menor de colunas para incluir no modelo.
#COLSET_01 = c("Plural", "Sex", "MomAge", "Weeks",
# "Gained", "Smoke",
# "BirthWeightOz")
# fit_weight_2 é o fModelo 4 ====
fit_weight_2 = lm(
BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke,
data = ncb_weight
)
summary(fit_weight_2)##
## Call:
## lm(formula = BirthWeightOz ~ Plural + Sex + MomAge + Weeks +
## Gained + Smoke, data = ncb_weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -63.672 -10.236 -0.157 10.486 49.522
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -63.14127 7.06144 -8.942 < 2e-16 ***
## PluralTwins -25.54412 2.71720 -9.401 < 2e-16 ***
## PluralTriplets -33.92719 8.45212 -4.014 6.28e-05 ***
## SexF -3.42455 0.88540 -3.868 0.000115 ***
## MomAge 0.50662 0.07338 6.904 7.65e-12 ***
## Weeks 4.16276 0.17797 23.390 < 2e-16 ***
## Gained 0.28536 0.03215 8.875 < 2e-16 ***
## SmokeY -7.24103 1.26076 -5.743 1.14e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.59 on 1401 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.4408, Adjusted R-squared: 0.438
## F-statistic: 157.8 on 7 and 1401 DF, p-value: < 2.2e-16
fit_weight_2b = lm(
BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke,
data = ncb_weight_NF)
summary(fit_weight_2b)##
## Call:
## lm(formula = BirthWeightOz ~ Plural + Sex + MomAge + Weeks +
## Gained + Smoke, data = ncb_weight_NF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -63.654 -10.185 -0.172 10.478 49.537
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -63.40360 7.05644 -8.985 < 2e-16 ***
## PluralTwins -25.56954 2.71478 -9.419 < 2e-16 ***
## PluralTriplets -33.98035 8.44453 -4.024 6.03e-05 ***
## SexF -3.38139 0.88489 -3.821 0.000139 ***
## MomAge 0.51259 0.07339 6.985 4.39e-12 ***
## Weeks 4.16531 0.17782 23.425 < 2e-16 ***
## Gained 0.28563 0.03213 8.891 < 2e-16 ***
## SmokeY -7.25583 1.25964 -5.760 1.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.57 on 1400 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.4419, Adjusted R-squared: 0.4391
## F-statistic: 158.3 on 7 and 1400 DF, p-value: < 2.2e-16
anova(fit_weight_1, fit_weight_2)## Analysis of Variance Table
##
## Model 1: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Marital + Gained +
## Smoke + Premie + Etnicidade
## Model 2: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1386 372384
## 2 1401 385454 -15 -13069 3.2429 2.512e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fit_weight_1b, fit_weight_2b)## Analysis of Variance Table
##
## Model 1: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Marital + Gained +
## Smoke + Premie + Etnicidade
## Model 2: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1386 372384
## 2 1400 384483 -14 -12098 3.2164 4.908e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Inesperadamente, o a soma dos quadrados dos resíduos aumentou quando as outras variáveis foram retiradas.
# fit_weight_5 é o fModelo 5 ====
if (useEtn) {
fit_weight_5 = lm(
BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke + Etnicidade,
data=ncb_weight)
sf5 = summary(fit_weight_5)
print(sf5)
anova(fit_weight_1, fit_weight_2, fit_weight_5)
} else {
fit_weight_5 = lm(
BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke + RaceMom,
data = ncb_weight)
sf5 = summary(fit_weight_5)
print(sf5)
# fit_weight_5{b,c} são variações do fModelo 3 ====
fit_weight_5b = lm(
BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke + MomRace,
data = ncb_weight)
summary(fit_weight_5b)
fit_weight_5c = lm(
BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke + HispMom,
data = ncb_weight)
summary(fit_weight_5c)
anova(fit_weight_1, fit_weight_2, fit_weight_5, fit_weight_5b, fit_weight_5c)
}##
## Call:
## lm(formula = BirthWeightOz ~ Plural + Sex + MomAge + Weeks +
## Gained + Smoke + Etnicidade, data = ncb_weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.805 -10.247 -0.376 10.312 52.716
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -61.07108 7.88766 -7.743 1.87e-14 ***
## PluralTwins -26.01446 2.70672 -9.611 < 2e-16 ***
## PluralTriplets -35.03222 8.40840 -4.166 3.29e-05 ***
## SexF -3.33780 0.88213 -3.784 0.000161 ***
## MomAge 0.44590 0.07632 5.842 6.40e-09 ***
## Weeks 4.12488 0.17780 23.199 < 2e-16 ***
## Gained 0.27935 0.03246 8.606 < 2e-16 ***
## SmokeY -7.66939 1.28020 -5.991 2.66e-09 ***
## Etnicidadeblack Black N -2.18985 3.72465 -0.588 0.556672
## Etnicidadeblack portoriq 11.56560 16.92197 0.683 0.494426
## Etnicidadeblack south american 5.59946 16.90440 0.331 0.740511
## Etnicidadecentro-south american 0.66689 5.05104 0.132 0.894979
## Etnicidadechinese 5.60070 12.22582 0.458 0.646949
## Etnicidadecuban 4.74089 12.21830 0.388 0.698064
## Etnicidadefilipino -29.30335 16.91408 -1.732 0.083410 .
## Etnicidadehispanic other 5.84796 10.20114 0.573 0.566558
## Etnicidademexican 2.61558 3.93672 0.664 0.506541
## EtnicidadeOtherAsianOrPacific -1.46203 5.05622 -0.289 0.772506
## Etnicidadeportoriq -12.73123 7.20327 -1.767 0.077377 .
## Etnicidadesouth americanIndian 9.05507 16.91154 0.535 0.592433
## Etnicidadewhite 2.50894 3.65806 0.686 0.492912
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.49 on 1388 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.4524, Adjusted R-squared: 0.4445
## F-statistic: 57.33 on 20 and 1388 DF, p-value: < 2.2e-16
## Analysis of Variance Table
##
## Model 1: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Marital + Gained +
## Smoke + Premie + Etnicidade
## Model 2: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke
## Model 3: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke +
## Etnicidade
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1386 372384
## 2 1401 385454 -15 -13069.3 3.2429 2.512e-05 ***
## 3 1388 377445 13 8008.8 2.2930 0.00536 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comparaCoeficientes(summary(fit_weight_5), sf1)A inclusão dos campos étnicos prejudicou o modelo, pois os erros cresceram ou ficaram confusos (p>0.05).
Houve diferença significativa entre os coeficientes do intercepto e do campo Weeks (mas não sei por que)
O modelo 5 chamou atenção para o campo Weeks. Pelo contexto, espera-se uma forte correlação entre o peso do nascituro e a idade gestacional do parto: quanto mais próximo do termo, maior o peso (o que acontece depois do termo?).
Esta seção tenta isolar a influência desse campo no modelo.
# fit_weight_W é um modelo linear simples para comparar com o fModelo 6 ====
fit_weight_W = lm(BirthWeightOz ~ Weeks, data=ncb_weight)
sfW = summary(fit_weight_W)
print(sfW)##
## Call:
## lm(formula = BirthWeightOz ~ Weeks, data = ncb_weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.480 -11.994 -0.286 11.908 58.006
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -73.1171 6.7817 -10.78 <2e-16 ***
## Weeks 4.9028 0.1752 27.99 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.99 on 1447 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.3512, Adjusted R-squared: 0.3508
## F-statistic: 783.4 on 1 and 1447 DF, p-value: < 2.2e-16
Um modelo que tem apenas o campo Weeks como preditor de BirthWeighOz explica pouco mais de 35% da variância. É bastante, comparado com o modelo 3, com todas as variáveis, que explica 45.9741839%
# anova(fit_weight_1, fit_weight_W, fit_weight_2)
print("Cadê o noNA?")## [1] "Cadê o noNA?"
if (useEtn) { # Campo "etnicidade" unificado no modelo 6 ====
fit_weight_6 = lm(BirthWeightOz ~ Plural + Sex + MomAge + Gained + Smoke + Etnicidade, data=ncb_weight)
sf6 = summary(fit_weight_6)
print(sf6)
anova(fit_weight_1, fit_weight_6) # compara variância de "1" com "6"
} else { # Campos étnicos separados (como no original) no modelo 6 ====
fit_weight_6 = lm(BirthWeightOz ~ Plural + Sex + MomAge + Gained + Smoke + RaceMom, data = ncb_weight)
sf6 = summary(fit_weight_6)
print(sf6)
fit_weight_6b = lm(BirthWeightOz ~ Plural + Sex + MomAge + Gained + Smoke + MomRace, data = ncb_weight)
sf6b = summary(fit_weight_6b)
print(sf6b)
fit_weight_6c = lm(BirthWeightOz ~ Plural + Sex + MomAge + Gained + Smoke + HispMom, data = ncb_weight)
sf6c = summary(fit_weight_6c)
print(sf6c)
anova(fit_weight_1, fit_weight_2, fit_weight_6, fit_weight_6b, fit_weight_6c)
}##
## Call:
## lm(formula = BirthWeightOz ~ Plural + Sex + MomAge + Gained +
## Smoke + Etnicidade, data = ncb_weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -94.525 -10.042 0.959 11.864 52.319
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92.92254 5.01749 18.520 < 2e-16 ***
## PluralTwins -45.55676 3.02917 -15.039 < 2e-16 ***
## PluralTriplets -67.39100 9.76461 -6.902 7.80e-12 ***
## SexF -2.27666 1.03741 -2.195 0.0284 *
## MomAge 0.54296 0.08974 6.050 1.86e-09 ***
## Gained 0.35434 0.03804 9.316 < 2e-16 ***
## SmokeY -8.92184 1.50624 -5.923 3.97e-09 ***
## Etnicidadeblack Black N -2.51296 4.38614 -0.573 0.5668
## Etnicidadeblack portoriq 8.28203 19.92673 0.416 0.6777
## Etnicidadeblack south american 5.05169 19.90672 0.254 0.7997
## Etnicidadecentro-south american -0.57218 5.94780 -0.096 0.9234
## Etnicidadechinese 1.43880 14.39566 0.100 0.9204
## Etnicidadecuban 6.84325 14.38796 0.476 0.6344
## Etnicidadefilipino -25.42861 19.91717 -1.277 0.2019
## Etnicidadehispanic other 4.25154 12.01266 0.354 0.7235
## Etnicidademexican 5.76263 4.63316 1.244 0.2138
## EtnicidadeOtherAsianOrPacific 0.72237 5.95321 0.121 0.9034
## Etnicidadeportoriq -8.50094 8.47990 -1.002 0.3163
## Etnicidadesouth americanIndian 12.94574 19.91416 0.650 0.5158
## Etnicidadewhite 3.57090 4.30742 0.829 0.4072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.42 on 1389 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.2401, Adjusted R-squared: 0.2297
## F-statistic: 23.09 on 19 and 1389 DF, p-value: < 2.2e-16
## Analysis of Variance Table
##
## Model 1: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Marital + Gained +
## Smoke + Premie + Etnicidade
## Model 2: BirthWeightOz ~ Plural + Sex + MomAge + Gained + Smoke + Etnicidade
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1386 372384
## 2 1389 523802 -3 -151417 187.86 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comparaCoeficientes(summary(fit_weight_6), sf1)Uma informação ausente sobre o peso do nascituro é: de qual bebê é o peso informado em BirthWeightOz? É a média dos dois ou três? É o menor (ou maior) deles? Essa falta de informação pode estar introduzindo alguma distorção no modelo. Façamos um sem gêmeos.
# fit_weight_7 é o fModelo 7 (sem gêmeos) ========
if (useEtn) { # Campo "etnicidade" unificado no modelo 7
fit_weight_7 = lm(
BirthWeightOz ~ Sex + MomAge + Gained + Smoke + Etnicidade,
data=ncb_weight)
sf7 = summary(fit_weight_7)
print(sf7)
anova(fit_weight_1, fit_weight_7, fit_weight_2) # compara variância de "1" com "2" e "7" =====
} else { # Campos étnicos separados (como no original) no modelo 7 =======
fit_weight_7 = lm(BirthWeightOz ~ Sex + MomAge + Gained + Smoke + RaceMom, data = ncb_weight)
sf7 = summary(fit_weight_7)
print(sf7)
fit_weight_7b = lm(BirthWeightOz ~ Sex + MomAge + Gained + Smoke + MomRace, data = ncb_weight)
sf7b = summary(fit_weight_7b)
print(sf7b)
fit_weight_7c = lm(BirthWeightOz ~ Sex + MomAge + Gained + Smoke + HispMom, data = ncb_weight)
sf7c = summary(fit_weight_7c)
print(sf7c)
anova(fit_weight_1, fit_weight_2, fit_weight_7, fit_weight_7b, fit_weight_7c)
}##
## Call:
## lm(formula = BirthWeightOz ~ Sex + MomAge + Gained + Smoke +
## Etnicidade, data = ncb_weight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -99.731 -9.860 1.728 12.816 54.469
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.47907 5.47331 17.627 < 2e-16 ***
## SexF -2.13674 1.13286 -1.886 0.0595 .
## MomAge 0.45340 0.09778 4.637 3.87e-06 ***
## Gained 0.29870 0.04136 7.222 8.43e-13 ***
## SmokeY -7.56524 1.64242 -4.606 4.48e-06 ***
## Etnicidadeblack Black N -3.62693 4.78926 -0.757 0.4490
## Etnicidadeblack portoriq 6.08428 21.76065 0.280 0.7798
## Etnicidadeblack south american 4.24173 21.73921 0.195 0.8453
## Etnicidadecentro-south american 0.18372 6.49515 0.028 0.9774
## Etnicidadechinese 3.09240 15.72049 0.197 0.8441
## Etnicidadecuban 6.84812 15.71250 0.436 0.6630
## Etnicidadefilipino -23.97521 21.75047 -1.102 0.2705
## Etnicidadehispanic other 4.87315 13.11846 0.371 0.7103
## Etnicidademexican 4.93014 5.05936 0.974 0.3300
## EtnicidadeOtherAsianOrPacific 1.48282 6.50105 0.228 0.8196
## Etnicidadeportoriq -8.37477 9.26055 -0.904 0.3660
## Etnicidadesouth americanIndian 15.19172 21.74690 0.699 0.4849
## Etnicidadewhite 1.96807 4.70272 0.418 0.6756
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.21 on 1391 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.0924, Adjusted R-squared: 0.08131
## F-statistic: 8.33 on 17 and 1391 DF, p-value: < 2.2e-16
## Analysis of Variance Table
##
## Model 1: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Marital + Gained +
## Smoke + Premie + Etnicidade
## Model 2: BirthWeightOz ~ Sex + MomAge + Gained + Smoke + Etnicidade
## Model 3: BirthWeightOz ~ Plural + Sex + MomAge + Weeks + Gained + Smoke
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 1386 372384
## 2 1391 625582 -5 -253197 188.48 < 2.2e-16 ***
## 3 1401 385454 -10 240128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comparaCoeficientes(summary(fit_weight_7), sf1) # , sf2, sf6)c. Verifique as premissas do modelo linear.
Os resíduos dos modelos 6 e 7 tem sua estatística W abaixo de 96%. Todos os outros ficara acima de 99%. Seus gráficos Q-Q apresentam deformidades que levantam suspeitas de que não sejam ~N.
source("verifica_residuos.R")
pro_res_1 = verifica_props_residuos(
fit_gained1$residuals,
ncb_births_noGainedNA$Gained,
fit_gained1$fitted.values
)
if (pro_res_1$is.small) {
print("Propriedades do modelo 1 OK")
} else {
print("Alguma propriedade do modelo 1 não tem o valor esperado:")
print_props_residuos(pro_res_1)
}## [1] "Propriedades do modelo 1 OK"
# bloco de código - item c
plot(predict(fit_gained1), resid(fit_gained1), # data=nc_births,
main="Valores ajustados versus resíduos para modelo Naïve")par(mfrow=c(2,2))
plot(fit_gained1)## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
par(mfrow=c(1,1))pro_res_2 = verifica_props_residuos(
fit_gained2$residuals,
ncb_births_noGainedNA$Gained,
fit_gained2$fitted.values
)
if (pro_res_2$is.small) {
print("Propriedades do modelo 2 OK")
} else {
print("Alguma propriedade do modelo 2 não tem o valor esperado:")
print_props_residuos(pro_res_2)
}## [1] "Propriedades do modelo 2 OK"
plot(predict(fit_gained2), resid(fit_gained2), # data=nc_births,
main="Valores ajustados versus resíduos para modelo 2 (aprimorado)")par(mfrow=c(2,2))
plot(fit_gained2)par(mfrow=c(1,1))TODO: Interpretação
par(mfrow = c(2, 2))
plot(fit_weight_1)## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
par(mfrow = c(2, 2))
plot(fit_weight_1b)## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383
## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383
O gráfico de resíduos x valores ajustados apresenta uma forte heterodascidade, com uma grande concentração em volta de 123 Oz.
par(mfrow = c(2, 2))
plot(fit_weight_5)## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
Weeks
par(mfrow = c(2, 2))
plot(fit_weight_6)## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
Plural
par(mfrow = c(2, 2))
plot(fit_weight_7)## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
## Warning: not plotting observations with leverage one:
## 1200, 1306, 1383, 1387
O gráfico de resíduos x valores ajustados apresenta uma forte heterodascidade, com uma grande concentração em volta de 123 Oz.
NADA POR AQUI
d. Proposta de modelo
Com base nas análises, proponha um ou mais modelos lineares multivariados. Explique a sua escolha.
# bloco de código - item d# bloco de código - item e
weeks = c(10, # Poucas chances de sobrevivência
11,
23, # 17% de sobrevivência
25,
39, # termo
39.5,
40.2,
44,
45,
46, # risco
47)
gained = c(0.5, 11.3, 43.8, 100, 110)
age = c(4, 6, 8, 10, 12, 14, # Abaixo
21.5, 29.3, # no meio
44, 45, 46, 60) # acima
underage_premie = list(Plural="Single", MomAge=10, Sex="M", Weeks=20, Gained=5, Smoke="N")
underage_tardie = list(Plural="Single", MomAge=10, Sex="M", Weeks=45, Gained=5, Smoke="N")
up1 = predict(fit_weight_2, underage_premie)
up2 = predict(fit_weight_2, underage_tardie)Quantas onças pesaria um bebê nascido nas circunstâncias:
| Modelo | underage/premie | underage/overtime | overage/premie | overage/tardie |
|---|---|---|---|---|
| 1 | 26.6069872 | 130.6760226 |
Em retrospecto, penso que teria perdido menos tempo se tivesse feito uma simples análise do R² ajustado de todas as variáveis numéricas (já que são poucas). Minha primeira escolha foi justamente a que tem menor potencial de explicação, com os dados disponíveis.
Esses valores sugerem um limite superior para o poder de explicação de um modelo linear que explique as obervações de cara uma dessas variáveis.
Os dados da tabela abaixo foram obtidos da execução do script varios_modelos.R, ideia que só tive perto do final do prazo para entregar o exercício.
| Coluna | R² ajustado |
|---|---|
| MomAge | 0.239690016396153 |
| Weeks | 0.588115108607321 |
| Gained | 0.0818566109518561 |
| BirthWeightOz | 0.568165398604105 |
Até agora, parece não haver evidências suficientes para compor um modelo linear do peso ganho em função de algumas das outras variáveis. O máximo possível de explicação é de 8% com todos os campos, e vários deles geram coeficientes sem significância estatística (\(p>5\%\)).
Quando passei a considerar o peso do bebê ao nascer (BirthWeightOz) como variável dependente, obtive modelos com maior poder de predição ($ {R_{ajust}}^2 > 40%$)